Introduction

Cashmere goat is a species known for its cashmere and has a stable heritability and a strong adaptability (Hu et al. 2012). Cashmere is a thin layer of fine wool that grows on the skin of goats and covers their coarse wool. It grows in winter and provides resistance against the wind and cold. In the spring, it begins to fall off which can adapt to climate change naturally. Cashmere is rare, special and most produced in the world, in which the yield of cashmere is generally high (Jin et al. 2018). Cashmere has great economic value, especially as a clothing product for warmth and comfort. With the popularity of the mass of cashmere, people's demand for cashmere is getting higher and higher. Although the output of cashmere is high, the overall trend of cashmere is too coarse, affecting the formation of textile crafts. The diameter of cashmere fiber plays a decisive role, which is the focus of our research.

In the past few years, the development of high-throughput analysis has led to widespread attention to miRNA because they are new regulators of cell development. MicroRNA is a non-coding small RNA with a length of about 18 to 25 nt, which acts as a negative regulator of mRNA expression by bundled to complementary sequences in the 3'-UTR of target mRNA and cause for translation inhibition or degradation of target (Lee et al. 1993). MiRNA is highly conservative among species. They are expressed in particular tissues at specific phases of development (Bashirullah et al. 2003). There is no doubt that miRNAs play a crux role in the process of wide range biological functions, such as development, differentiation, tissue morphogenesis, and even some diseases (Ma et al. 2018). A series of studies have shown that many miRNA may be involved in the growth of cashmere goats’ skin (Liu et al. 2012). For example, (Zhang et al. 2007) demonstrated that mmu-miR-720 and mmu-miR-199b were expressed primarily in the skin tissue of cashmere goats. Bai et al. (2016) detected that the expressions of miR-103-3p, miR-15b-5p, miR-17-5p, miR-30c-5p, miR-200b, miR-199a-3p, miR-199a-5p, miR-30a-5p, and miR-29a-3p were significantly different between anagen and telogen skin in LCG. Liu et al. (2012) confirmed that the veridicality of 316 known miRNAs and 22 novel miRNAs in goat skin. On the other hand, the expression characteristics in cashmere skin have not been fully understood, so they can still be used for further investigation.

In our study of the current, to find miRNA that has a potential effect on cashmere fineness in Cashmere goat skin, we use a technology of high-throughput sequencing to judge the miRNAs in LCG goat skin sample, and numerous of miRNAs were acquired in goat skins. To further explore the potential role of miRNAs in cashmere fineness, we build the regulatory network of miRNAs-miRNAs, and Protein-Protein Interaction (PPI) in LCG skin. Thus, this new miRNA study will lay a firm base to further explore its role in regulating skin.

Materials and Methods

Ethics

All animal studies have been conducted and approved by the Shenyang Agricultural University Animal Laboratory Committee (No. 201806011, Shenyang, China).

Goat skin samples

Skin samples of three female LCGs were collected. All the source animals were maintained under the same conditions. In handling the animals, we attempted to minimize stress and used local anesthesia (procaine) to avoid pain and reduce blood circulation. All experimental steps in this study have been approved and implemented based on the Experimental Animal Management Committee of Shenyang Agricultural University (201306007). For the analysis of miRNAs in the skin, we selected one-year-old goats, which are in the anagen of the October period of hair growth, for skin sampling. The skin samples were collected from the shoulders and hip, where cashmere grows densely in Cashmere goats and sealed in cryogenic tubes containing TRIzol reagent. The test tubes were immediately stored in liquid nitrogen for subsequent RNA isolation.

RNA isolation, library preparation, and sequencing

RNA degradation and contamination were evaluated on a 1% agarose gels. Is RNA purity by using NanoPhotometer® Spectrophotometer (IMPLEN, U.S.A.). RNA concentration was measured by using the Qubit® RNA Assay Kit with the Qubit® 2.0 Fluorometer (Life Technologies, C.A., U.S.A.). RNA integrity was assessed using the RNA Nano 6000 assay kit (Agilent Technologies, CA, USA) on the Bioanalyzer 2100 system. After the samples were tested, the library was built using the miRNA sample preparation kit. Raw image data files obtained using the Illumina HiSeq^TM 2000/MiSeq sequencing platform were converted to the original sequencing sequences (sequenced reads) and screened for specific lengths of sRNA for subsequent analysis. We filtered the raw reads to exclude low-quality sequences, 5' joint contamination reads, and poly-A/T/G/C reads. Then, Bowtie software was used to map the length-selected sRNAs to the reference genome and analyse the distribution of small RNAs on the genome. The sRNA tags were mapped to the RepeatMasker and Rfam databases to remove tags from protein-coding genes, repetitive sequences, rRNA, tRNA, snRNA, and snoRNA.

Alignment of known miRNAs and prediction of novel miRNAs

MiRNAs have been identified to abide by the work process. First, Tophat2 was used to map reads to a reference genome. Assemble the mapped readings of each sample by cufflinks. T Thereafter, transcripts that did not meet the criteria were rejected (length ≥ 200 nt, exon number ≥ 2 and reading coverage ≥ 3). In the case of PFAM, each transcript is translated into all three possible reading frames and the presence of known protein family domains recorded in the Pfam database (http://pfam.xfam.org) is determined using a Pfam scan, was identified. Reference sequences were obtained using Mir Base 20.0, and potential miRNAs were obtained using the software of mirdeep2 (Friedlander et al. 2012) and the sRNA-tools-cli to draw secondary structures. The properties of the hairpin structure of the miRNA precursor can be used to predict novel miRNAs. The miREvo (Wen et al. 2010) and miRDeep2 software packages were integrated to predict new miRNAs by examining the secondary structure, the Dicer cleavage site and the minimum free energy of the small RNA tags that were not in the previous steps were annotated.

MiRNA family analysis

We investigated the phenomenon of miRNA families identified from samples from other species. In our analysis process, miFam.dat (http://www.mirbase.org/ftp.shtml) was used to examine known miRNAs to find families and submit new miRNA precursors to Rfam (http: // rfam.sanger). transmission. ac.uk/search/) search the Rfam family.

MiRNA target gene prediction

The prediction of the target genes of miRNAs was carried out by miRanda (Enright et al. 2003) or miRDB (Score> 85) (http://mirdb.org/) for animals. The combined results of the two programs were considered. R language program was used to analyze the difference between LCG and Inner Mongolia Cashmere goat (MCG) (Liu et al. 2012) data (q-value ≤ 0.05 and |log 2-Fold Change| ≥ 1).

GO and KEGG enrichment analyses

A GOseq-based, non-central and hypergeometric Wallenius distribution (Young et al. 2010) was implemented for the GO enrichment analysis, which can compensate for the gene length distortion. KEGG (Kanehisa et al. 2008) is a database reservoir for understanding the features and functions of a biological system (such as cells, organisms, or ecosystems) based on molecular-level information, which is based on sequencing of molecular data sets generated by the genome and other experimental high-throughput technologies (http://www.genome.jp/kegg/). We used the software KOBAS (Smoot et al. 2011) to statistically test the enrichment of the KEGG target gene pathway.

PPI network construction

The interactions between the target genes generated in the miRNA analysis were identified in the STRING database, and the interaction data were imported into Cytoscape 3.5.1 (Mao et al. 2005) software to visualize the interaction network. Genes with a degree > 20 were considered as hub genes (proteins).

RNA preparation and quantitative real-time PCR

The previous results found that the SNP of the NFKBIA gene was significantly correlated with cashmere fineness, and the eexperimental LCG population was divided into three genotype populations with different fineness according to the SNP locus of NFKBIA gene, and then qPCR was verified with MCG. Total RNA was extracted from 9 LCGs and 3 MCGs goat skin samples, assessed its quality and concentration. Samples that showed reasonable RNA quality were selected for further analysis. Twenty-seven miRNAs were selected at random and verified with qPCR. The primer sequences were designed using Primer Premier 5.0. MiR-let-7a was used as the standard to verify the results. At least 3 samples were used in each analysis, and each analysis was repeated three times. The 2^−ΔΔCt method was used to analyse the relative expression levels of qPCR. The miRNA levels, normalized against 18S, were evaluated by the 2^−ΔΔCt method.

Results

RNA extraction and sequencing

We examined the expression of miRNAs and small RNAs in LCG by collecting skin samples, extracting RNA and sequencing miRNAs. After sequencing, we obtained average 15245657 raw reads from the skin samples. These raw reads were aligned to the reference gene and annotated to confirm the quality of the samples. After filtered out the low-quality reads, 14801303 clean reads were obtained. It accounts for 97.09% of the total readings. The lengths of the sRNA range of between 18 nt and 35 nt. Reads of 22 nt were the most abundant of the clean reads, accounting for 42.95% of the clean reads. The second and third most common read lengths were 21 nt (17.22%) and 23 nt (13.62%), respectively. These results showed the reads were small RNAs (Fig. 1). Then, the reads were aligned to the reference genome, resulting in 6391961 mapped reads, which accounts for 47.10% of the total reads. The sample reads were then mapped to the same chain in the direction of the reference sequence, resulting in 4525684 mapped reads, accounting for 33.35% of the total reads. There were 447744 reads of unique mapped small RNAs, which were mapped to the opposite direction of the reference sequence; 1866277 reads were obtained, accounting for 13.75% of the total reads (Table 1). Density statistics relative to the chromosomes on the genome were calculated to determine the distribution of reads on each chromosome (Fig. 2).

Identification and expression analysis of miRNAs

To recognize known miRNAs in the samples, we compared the sample data to a range of sequences in miRBase to get sRNA details for each sample match (Table 2). The total number of sRNAs that matched all of the samples to known miRNA precursors was 4894375, and the number of species was 4911. The process of miRNA development from precursor to mature body involves digestion by Dicer; the specificity of the cleavage site results in a strong bias in the first base of the miRNA mature sequence. Therefore, the first base distribution of miRNAs of different lengths was analysed (Fig. 3A) in addition to the base distribution statistics of miRNAs (Fig. 3B). The expression level and mature sequences of the top 10 miRNAs are presented in Table 3. Among the known miRNAs, miR-26a exhibited the highest expression level.

The underlying principle of miRNA prediction involves the highly conserved module sequence upstream of the hairpin structure. MiRNA is located on the two arms of the hairpin structure. The restriction site of the enzyme cutting site is conserved, and the migration does not typically exceed 3 nt. The precursors of the hairpin have low free folding energy. Characteristic hairpin structures of miRNA precursors that can be used to predict novel miRNAs have been analysed. As a result, we predicted a total of 45 new miRNAs (Table 4). After the prediction of new miRNAs in the skin, the first base distribution of miRNAs of different lengths was obtained (Fig. 4A). Also, the base distribution statistics of each point of miRNA were obtained (Fig. 4B). The top five novel miRNAs Table 1: High-throughput data classification

Table 2: Details of each sample matching sRNA

Table 3: The 10 most abundantly expressed miRNAs in LCG

Table 4: Predicted novel miRNAs and sample sRNA comparison statistics

Table 5: The 5 most abundantly expressed novel miRNAs in LCG

are shown in Table 5, of which 27 novel miRNAs exhibited the highest level of expression.

MiRNA family analysis

Family analysis of the detected known and novel miRNAs was conducted, and the presence of the miRNA families in other species was investigated (Fig. 5). To explore the extent of this evolution, we performed comparisons with 206 species, including Bos taurus, Ovis aries, Sus scrofa, Equus caballus, Pan troglodytes and Mus musculus. The number of known and new miRNA precursors was 182. miR124 was found in 77 species, and miR10 was found in 76 species, indicating that miR124, miR9, and miR10 exist widely among species. The greatest miRNA family is the miR-2284 family, which has 71 family members, followed by the miR-154 family, which has 33 family members, and the miR-379 and miR-10 families, which have 12 family members. Other families had only one family member.

Fig. 1: The distribution of fragments of miRNAs of different lengths

Fig. 2: The distribution density map of reads on each chromosome. Note: The outermost circle in the figure represents the chromosomes selected for display; the grey background region in the middle is the distribution of 10,000 reads, with the red reads mapping to the positive chain and the blue reads to the negative chain. The innermost circle provides a comparison of all reads on the chromosome; orange denotes positive chain coverage, and green denotes negative chain coverages. Singular points that exceed the mean of all cover sets+3 standard deviations were excluded

GO and KEGG enrichment analyses

We explored the potential biological features of miRNAs in the skin of Cashmere goats by targeting their target genes. The GO enrichment analysis of all predicted target genes was performed and revealed that the target gene functions mainly involved GO term 0003463 "structural constituent of ribosome" and those associated with the activity of DNA polymerase, RNA binding, the activity of RNA-directed DNA polymerase, and molecular function.

By GO analysis, we identified 182 (P < 0.05) significant enrichments of these top gene signatures, which were classified into 3 GO categories: biological process (BP, 101), molecular function (MF, 61) and cellular component (CC, 20). The top 40 target genes that were significantly enriched in each GO were counted, and the outcomes were exhibited as the histogram (Fig. 6). Main pathway enrichment can define key biochemical metabolic pathways and signal transduction pathways involved in candidate target genes. The total number of KEGG pathways showing

Fig. 3: (A) First base preferences of known miRNAs 18-30 nt in length. (B) Base preferences of known miRNAs.

Fig. 5: Number of miRNAs of each family

Fig. 6: Candidate target gene GO enrichment histogram.

Fig. 4: (A) First base preferences of new miRNAs 18-30 nt in length. (B) Base preferences of novel miRNAs

enrichment was 267. KEGG enrichment is measured by rich factors, Q values and the number of genes enriched in the said pathway. The higher the Rich factor, the higher the degree of enrichment. Analysis of the KEGG pathways indicated that four major biological pathways (Endocytosis pathway, Progesterone-mediated oocyte maturation, Wnt signalling pathway, and PI3K-Akt signalling pathway) may play an important role in hair circulation (Fig. 7).

MiRNA Target gene prediction

Fig. 7: Candidate target gene KEGG enrichment scatter plot

Fig. 8: Network Diagram of 27 miRNAs and their Target genes. Green: miRNAs, Blue: Target genes.

To predicted target genes because of the known and novel miRNAs after the determination of the correlation between miRNAs and target genes. The number of known and novel miRNAs was 509, and a total of 21361 target genes were predicted. We compared the target genes of the 27 miRNAs with our data using the online software miRDB (Fig. 8) and selected the three target genes with the highest degrees of gravity: CHRM2, FGL2, and GPR6. MiRNA mainly acts on target genes at the post-transcriptional level, and all three genes act on related regulatory factors and signalling pathways, thus affecting the growth and development of skin. We selected one of the pathways for the next target gene based on the functions provided by KEGG. Through the comparison of our high-throughput data with the MCG data (Liu et al. 2012) and the involvement of candidate miRNAs in the endocytosis pathway, we found 237 significant expressed miRNAs,147 differentially expressed miRNAs comprising 69 upregulated and 78 downregulated miRNAs (q-value ≤ 0.05 and |log 2 Fold Change| ≥ 1), and their related target genes were obtained (Fig. 9). In addition, the 27 miRNAs shown in Figure 8 were found to interact with target ribonucleic acids and were subjected to subsequent qPCR analysis. The corresponding miRNAs target gene will be the main candidate for further study of cashmere fineness.

Fig. 9: The network of differentially expressed miRNAs. Green: Downregulated miRNAs, Red: Upregulated miRNAs, Blue: Target genes.

Fig. 10: Protein-protein interaction network. Dashed lines represent the shortest paths, and the associated label denotes the length. The darker the node colour, the higher the score

Protein-protein interaction analysis

For species with annotation information in the protein-protein interaction (PPI) database, we performed a target gene protein interaction network analysis of miRNAs. For the list of target genes generated in the miRNA analysis, we identified the interaction relationship between these target genes in the STRING database (string score≥ 700). Using Cytoscape 3.5.1, an interaction network with 2040 nodes and 3645 edges was obtained. A gene with a degree larger than the threshold value (degree = 20) was considered a hub gene. Twenty-eight genes in the PPI network were identified as hub genes (Fig. 10); these might play important roles in the biological processes influencing cashmere fineness. The gene PCNA showed the highest degree (degree = 47) in the network, followed by RPL11 (degree = 40) and RPS13 (degree = 33).

Validation of miRNAs by qPCR

To affirm the RNA-seq and expression levels, we performed further analyses of LCG and MCG skin (Liu et al. 2012) tissue. Twenty-seven miRNAs co-existing in LCG and MCG skin (Liu et al. 2012) were randomly selected from different genotypes for validation (Fig. 11). The qPCR results showed that 27 miRNAs and LCG_GG types negatively regulate the fineness of cashmere. Among them, miR-101 has the least negative regulation, and miR-

Fig. 11: qPCR validation of miRNAs. Red: high expression, Blue: low expression. The degree of colour saturation reflects the level of expression

30e-5p and miR-30a-5p have the most negative regulation of cashmere fineness. The relative expression levels of 27 different miRNAs were compatible with the results of the solexa sequence, exhibiting similar trends, which confirmed that these miRNAs exist in LCG skin and may impact cashmere fineness. These results revealed that these miRNAs might be critical in Cashmere goat skin.

Discussion

The Cashmere goat is a species with high cashmere yield and high economic value. The skin of this goat is the key factor determining the cashmere yield. As one of the largest populations of sheep, the Cashmere goat population has made a huge contribution to China (Su et al. 2015). MiRNA is a non-coding RNA (ncRNA) that regulates gene expression (Bernardo et al. 2015). To date, approximately 25000 miRNAs have been confirmed in 193 animals, plant and microorganism species (Alvarezgarcia and Miska 2005). Over the past few years, it has been increasingly shown to play a vital role in the miRNAs is in various tissues, such as the brain (Mohammed et al. 2017), heart (Bernardo et al. 2015), liver (Roderburg and Trautwein 2016) and uterus (Ahn et al. 2017). However, very few researches focused on the hidden part of the miRNA in the skin’s growth of goat.

In the current study, miRNAs and mRNA of the skin of the Cashmere goat were systematically associated, using RNA-sequencing technology. We subsequently characterized the putative miRNAs to elucidate their various characteristics in order to provide an overview of the uncharacterized and the relationship with the development of goat skin. In the skin of Cashmere goat, a mountain of non-coding RNAs has been discovered (Liu et al. 2012; Schmitt et al. 2012; Yuan et al. 2013). Among these miRNAs, we identified a total of 464 known miRNAs and 45novel miRNAs from our high-throughput sequencing results.

Family analyses based on high-throughput sequencing identified 182 families. In our results, miR-127, miR-154, miR-229, miR-379, miR-368, and miR-432 are all conserved in goat skin. Among these miRNAs have also been shown to be conserved in mammals, in which, miR-368 family is only expressed in sheep and goats, not in cattle, horses, and pigs (Liu et al. 2012). MiR-368 may be associated with skin development and normal apoptosis, so these miRNAs may be involved in regulating goat skin development.

We then analyzed relevant gene functions and their corresponding signaling pathways through enrichment analysis of GO and KEGG terms. A number of the related candidate target genes were related to the molecular function term (2588), accounting for 89.8% of the total number of candidate target genes on the annotation. Our results suggest that several pathways may affect the regulation of cashmere. There is increasing evidence that hundreds of miRNAs affect most transcriptomes and a wide range of biological processes and signalling pathways, such as the Wnt, Notch and bone morphogenetic protein (BMP) signalling pathways (Assa-Kunik et al. 2007; Mardaryev et al. 2010; Bai et al. 2018). It has been reported that the development of wool follicles is affected by the WNT pathway (Park et al. 2012). Other enriched terms included the metabolic pathways in cancer, the MAPK signalling pathway, endocytosis, and focal adhesions, indicating that these pathways are likely to play an important role during the hair cycle (Yuan et al. 2013). The PI3K-AKT pathway has been shown to be involved in the treatment of skin cancer (Macdonald et al. 2015) and is critical for cellular growth and metabolism, including hair growth (Suzuki et al. 2017). Our data also enriched these pathways; it further illustrates the importance of these miRNAs in goat skin.

A total of 21361 miRNA genes have been identified in LCG skin tissue. Some genes, such as HOXC850, BMP451, SDHA, YWHAZ, and UBC52, have been reported to be associated with the length and fiber in LCG. Ten genes are associated with cashmere fineness: AKT1, ALX4, HK1, NT-3, POLD2, PSMA2, MAPK8, RYR3, VPS39, and STARD9. GWAs have suggested that these genes are associated with cashmere fineness in Inner Mongolia Cashmere goats (Qiao et al. 2017). HF genes linked to the cycle, some like MSX2, WNT3A, and HOXC13, have been discovered (Wang et al. 2017). Thus, the target gene is our main research direction in the future.

Further investigation of miRNA gene interactions, we use the miRDB tool to obtain some target genes that may have potential effects. The results indicated that three targets are higher enrichment; we surmise that they may have some latent effects. From a functional point of view, a single miRNA can regulate the expression of more than 100 target mRNAs, and vice versa, the expression of the mRNA can be regulated by many different miRNAs (Lewis et al. 2003; Kim 2005). Three target genes and their associated miRNAs can be clearly seen from the network. As an example, the three-target gene (CHRM2, FGL2, and GPR6) was all involved in miR-93, in which CHRM2, FGL2 all engaged in miR-17-5p, and FGL2, GPR6 also shared in miRNA-106b. There was evidence that miR-93 and miR-17-5P were upregulated in the catagen than anagen stage of skin in Cashmere goat (Yuan et al. 2013). Therefore, it proved potential reference value for cashmere fiber fineness and the expression analysis of miRNAs in LCG.

PPIs are important components of the whole cell-cell interaction system in living cells. The PPI is considered a functional unit of most biological processes, including the development and metastasis of human cancer (Zhang et al. 2017). PPI network analysis showed that PCNA, RPL11, and RPS13 had high degrees in the network. Proliferating cell nuclear antigen (PCNA) was previously thought to be a DNA sliding fixture for DNA polymerase replication and is an important component of eukaryotic chromosomal DNA replicons. However, subsequent studies have shown that it has significant capabilities to interact with multiple partners involved in multiple metabolic pathways, including Okazaki fragment processing, DNA repair, trans-damage DNA synthesis, DNA methylation, chromatin remodeling, and the cell cycle regulation (Maga and Hubscher 2003). Polymorphisms and their possible associations with cashmere fiber traits in the 5' upstream region of the prolactin receptor gene (PRLR) (5'UTR) have been reported, and PRLR encodes anterior pituitary hormones involved in various physiological activities (Zhou et al. 2011). RPL11 is involved in the downregulation of c-Myc, a transcription factor that plays a key role in chromosomal nucleolar biogenesis. Studies have shown that the binding of L11 to Mdm2 at the promoter results in Mdm2-mediated remission of p53 transcriptional repression (Mahata et al. 2012). Hybrid mutation of ribosomal protein L11 (RPL11) causes Diamond-Black fan anaemia-7 (DBA7) (Carlston et al. 2017). RPL11-deficient zebrafish contain a retroviral insertion in the first intron that interferes with mRNA splicing, resulting in reduced protein expression (Amsterdam et al. 2004).

The ribosomal protein S13 mRNA (RPS13) serves as a stable, moderately abundant reference gene for qPCR to normalize HIV/SIV RNA copy numbers (Robichaux et al. 2016). The RPS13 gene appears to play an important role in the development of NK/T-cell lymphoma (Yang et al. 2006). Studies have shown that the RPS13 gene encoding the ribosomal protein S13 and the ATP6 gene encoding the ATP synthase complex subunit 6 are co-transcribed (Takvorian et al. 1997). Therefore, we speculate that the study of PPI can be helpful for understanding the biological mechanisms related to the regulation of cashmere fineness and may be an effective method for predicting the genes related to the regulation of cashmere fineness.

To validate the sequencing data, we randomly selected 27 miRNAs and divided them into 3 different genotypes and verified expression levels in LCG by qPCR. The qPCR results of miRNAs are all negative regulation of cashmere fineness. Among them, miR-101, miR-30e-5p, and miR-30a-5p are thought to have potential effects on regulating skin. Some miRNAs have been reported to be expressed on Cashmere goats (Li et al. 2016). For instance, miR-30a-5p, miR-30e-5p significantly up-regulated on Cashmere goat skin during the catagen and/or telogen stages (Assa-Kunik et al. 2007). In our result, the differentially expressed miR-101 is involved in the endocytosis pathway, which is associated with cashmere fineness. Studies have shown that miR-101 is involved in a variety of cancer-related biological processes, including proliferation, apoptosis, angiogenesis, drug resistance, invasion and metastasis (Wang et al. 2018). Based on our results, it can be concluded that these three kinds of miRNAs (miR-101, miR-30e-5p, and miR-30a-5p) may be involved in the control of Cashmere goat skin.

Conclusion

This research offers valuable information for the role of miRNAs in regulating the development of cashmere fineness phenotype. 464 known miRNAs and 45 new miRNAs were found in the skin. The potential functions of the target genes were identified by KEGG enrichment assay and indicate the potential roles which miRNAs in skin, hair follicle development and growth. The expression level of miRNAs in the LCG_GG genotype is higher than that of the other genotypes, and the fineness of the cashmere can be inhibited. We hypothesize that miR-101, miR-30e-5p, and miR-30a-5p may have a regulatory effect on Cashmere goat skin. The results of this research are helpful to further elucidate the role of these miRNAs in the phenotypic differential expression of LCG skin.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (No. 31872325, 31802038, 31672388), the Liaoning Province Foundation of Natural Science project, China (No. 2015020758).

References

Ahn SH, V Singh, C Tayade (2017). Biomarkers in endometriosis: Challenges and opportunities. Fert Steril 107:523‒532

Alvarezgarcia I, EA Miska (2005). MicroRNA functions in animal development and human disease. Development 132:4653‒4662

Amsterdam A, RM Nissen, Z Sun, EC Swindell, S Farrington, N Hopkins (2004) Identification of 315 genes essential for early zebrafish development. Proc Natl Acad Sci USA 101:12792‒12797

Assa-Kunik E, IL Torres, ED Schejter, DS Johnston, BZ Shilom (2007). Drosophila follicle cells are patterned by multiple levels of notch signaling and antagonism between the notch and JAK/STAT pathways. Development 134:1161‒1169

Bai WL, SJ Zhao, ZY Wang, YB Zhu, YL Dang, YY Cong, HL Xue, W Wang, L Deng, D Guo, SQ Wang, YX Zhu, RH Yin (2018). LncRNAs in secondary hair follicle of Cashmere goat: Identification, expression, and their regulatory network in wnt signaling pathway. Anim Biotechnol 29:199‒211

Bai WL, YL Dang, RH Yin, WQ Jiang, ZY Wang, YB Zhu, SQ Wang, YY Zhao, L Deng, GB Luo, SH Yang (2016). Differential expression of micrornas and their regulatory networks in skin tissue of Liaoning Cashmere goat during hair follicle cycles. Anim Biotechnol 27:104‒112

Bashirullah A, AE Pasquinelli, AA Kiger, N Perrimon, G Ruvkun, CS Thummel (2003). Coordinate regulation of small temporal RNAs at the onset of Drosophila metamorphosis. Dev Biol 259:1‒8

Bernardo BC, JY Ooi, RC Lin, JR McMullen (2015). MiRNA therapeutics: A new class of drugs with potential therapeutic applications in the heart. Future Med Chem 7:1771‒1792

Carlston CM, ZA Afify, JC Palumbos, H Bagley, C Barbagelata, WL Wooderchak-Donahue, R Mao, JC Carey (2017). Variable expressivity and incomplete penetrance in a large family with non-lassical Diamond-Blackfan anemia associated with ribosomal protein L11 splicing variant. Amer J Med Genet A 173:2622‒2627

Enright AJ, B John, U Gaul, T Tuschl, C Sander, DS Marks (2003). MicroRNA targets in Drosophila. Genome Biol 5:1–14

Friedlander MR, SD Mackowiak, N Li, W Chen, N Rajewsky (2012). MiRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucl Acids Res 40:37‒52

Hu PF, WJ Guan, XC Li, WX Zhang, CL Li, YH Ma (2012). Study on characteristics of in vitro culture and intracellular transduction of exogenous proteins in fibroblast cell line of Liaoning Cashmere goat. Mol Biol Rep 40:327‒336

Jin M, JY Zhang, MX Chu, J Piao, JA Piao, FQ Zhao (2018). Cashmere growth control in Liaoning cashmere goat by ovarian carcinoma immunoreactive antigen-like protein 2 and decorin genes. Asian-Aust J Sci 31:650‒657

Kanehisa M, M Araki, S Goto, M Hattori, M Hirakawa, M Itoh, T Katayama, S Kawashima, S Okuda, T Tokimatsu, Y Yamanishi (2008). KEGG for linking genomes to life and the environment. Nucl Acids Res 36:480‒484

Kim VN (2005). Small RNAs: Classification, biogenesis, and function. Mol Cells 19:1‒15

Lee RC, RL Feinbaum, V Ambros (1993). The C elegansheterochronic gene lin-4 encodes small RNAs with antisense complementarity to lin-14. Cell 75:843‒854

Lewis BP, IH Shih, MW Jones-Rhoades, DP Bartel, CB Burge (2003). Prediction of mammalian microRNA targets. Cell 115:787‒798

Li J, H Qu, H Jiang, Z Zhao, Q Zhang (2016). Transcriptome-wide comparative analysis of microRNA profiles in the Telogen skins of Liaoning Cashmere goats (Capra hircus) and fine-wool sheep (Ovis aries) by solexa deep sequencing. DNA Cell Biol 39:696‒705

Liu Z, H Xiao, H Li, Y Zhao, S Lai, X Yu, T Cai, C Du, W Zhang, J Li (2012). Identification of conserved and novel microRNAs in cashmere goat skin by deep sequencing. PLoS One 7; Article e50001

Ma T, J Li, Q Jiang, S Wu, H Jiang, Q Zhang (2018). Differential expression of miR-let7a in hair follicle cycle of Liaoning cashmere goats and identification of its targets. Funct Integr Genomics 18:701‒707

Macdonald JB, B Macdonald, LE Golitz, P LoRusso, A Sekulic (2015). Cutaneous adverse effects of targeted therapies: Part II: Inhibitors of intracellular molecular signaling pathways. J Amer Acad Dermatol 72:221‒236

Maga G, U Hubscher (2003). Proliferating cell nuclear antigen (PCNA): A dancer with many partners. J Cell Sci 116:3051‒3060

Mahata B, A Sundqvist, DP Xirodimas (2012). Recruitment of RPL11 at promoter sites of p53-regulated genes upon nucleolar stress through NEDD8 and in an Mdm2-dependent manner. Oncogene 31:3060‒3071

Mao X, T Cai, JG Olyarchuk, L Wei (2005). Automated genome annotation and pathway identification using the KEGG orthology (KO) as a controlled vocabulary. Bioinformatics 21:3787‒3793

Mardaryev AN, MI Ahmed, NV Vlahov, MY Fessing, JH Gill, AA Sharov, NV Botchkareva (2010). Micro-RNA-31 controls hair cycle-associated changes in gene expression programs of the skin and hair follicle. J Feder Amer Soc Exp Biol 24:3869‒3881Mohammed CPD, JS Park, HG Nam, K Kim (2017). MicroRNAs in brain aging. Mech Ageing Dev 168:3‒9

Park PJ, BS Moon, SH Lee, SN Kim, AR Kim, HJ Kim, WS Park, KY Choi, EG Cho, TR Lee (2012). Hair growth-promoting effect of Aconiti ciliare tuber, extract mediated by the activation of wnt/β-catenin signaling. Life Sci 91:935‒943

Qiao X, R Su, Y Wang, R Wang, T Yang, X Li, W Chen, S He, Y Jiang, Q Xu, W Wan, Y Zhang, W Zhang, J Chen, B Liu, X Liu, Y Fan, D Chen, H Jiang, D Fang, Z Liu, X Wang, Y Zhang, D Mao, Z Wang, R Di, Q Zhao, T Zhong, H Yang, J Wang, W Wang, Y Dong, X Chen, X Xu, J Li (2017). Genome-wide target enrichment-aided chip design: A 66K SNP chip for cashmere goat. Sci Rep 7; Article 8621

Robichaux S, N Lacour, GJ Bagby, AM Amedee (2016). Validation of RPS13 as a reference gene for absolute quantification of SIV RNA in tissue of rhesus macaques. J Virol Meth 236:245‒251

Roderburg C, C Trautwein (2016). Cell-specific functions of miRNA in the liver. J Hepatol 66:655‒656

Schmitt MJ, D Philippidou, SE Reinsbach, C Margue, A Wienecke-Baldacchino, D Nashan, I Behrmann, S Kreis (2012). Interferon-g-induced activation of signal transducer and activator of transcription 1 (STAT1) up-regulates the tumor suppressing microRNA-29 family in melanoma cells. Cell Commun Signal 10:41‒55

Smoot ME, K Ono, J Ruscheinski, PL Wang, T Ideker (2011). Cytoscape 28 New features for data integration and network visualization. Bioinformatics 27:431‒432

Su R, S Fu, Y Zhang, R Wang, Y Zhou, J Li, W Zhang (2015). Comparative genomic approach reveals novel conserved microRNAs in Inner Mongolia cashmere goat skin and longissimus dorsi. Mol Biol Rep 42:989‒995

Suzuki Y, Y Enokido, K Yamada, M Inaba, K Kuwata, N Hanada, T Morishita, S Mizuno, N Wakamatsu (2017). The effect of rapamycin, NVP-BEZ235, aspirin, and metformin on PI3K/AKT/mTOR signaling pathway of PIK3CA-related overgrowth spectrum (PROS). Oncotarget 8:45470‒45483

Takvorian A, JL Coville, N Haouazine-Takvorian, A Rode, C Hartmann (1997). The wheat mitochondrial rps13 gene: RNA editing and co-transcription with the atp6 gene. Curr Genet 31:497‒502

Wang CZ, F Deng, H Li, DD Wang, W Zhang, L Ding, JH Tang (2018). MiR-101: A potential therapeutic target of cancers. Amer J Transl Res 10:3310‒3321

Wang S, W Ge, Z Luo, Y Guo, B Jiao, L Qu, Z Zhang, X Wang (2017). Integrated analysis of coding genes and non-coding RNAs during hair follicle cycle of cashmere goat (Capra hircus). BMC Genomics 18; Article 767

Wen M, Y Shen, S Shi, T Tang (2010). MiREvo: An Integrative MicroRNA Evolutionary Analysis Platform for Next-generation Sequencing Experiments. BMC Bioinform 13:140‒150

Yang F, WP Liu, MX He, QL Tang, S Zhao, WY Zhang, QJ Xia, GD Li (2006). Real-time fluorescence quantitative PCR in detecting ribosome protein S13 (RPS13) gene expression in NK/T cell lymphoma. J Sichuan Univ 37:464‒466

Young MD, MJ Wakefield, GK Smyth, A Oshlack (2010). Goseq: Gene ontology testing for RNA-seq accounting for selection bias. Genome Biol 11; Article R14

Yuan C, X Wang, R Geng, X He, L Qu, Y Chen (2013). Discovery of cashmere goat (Capra hircus) microRNAs in skin and hair follicles by Solexa sequencing. BMC Genomics 14; Article 511

Zhang J, TD Le, L Liu, J Li (2017). Inferring miRNA sponge co-regulation of protein-protein interactions in human breast cancer. BMC Bioinform 18; Article 243

Zhang WG, JH Wu, JQ Li, Y Midori (2007). A Subset of Skin-Expressed microRNAs with Possible Roles in Goat and Sheep Hair Growth Based on Expression Profiling of Mammalian microRNAs. OMICS J Integr Biol 11:385‒396

Zhou JP, XP Zhu, W Zhang, F Qin, SW Zhang, ZH Jia (2011). Short Communication A novel single-nucleotide polymorphism in the 5' upstream region of the prolactin receptor gene is associated with fiber traits in Liaoning cashmere goats. Genet Mol Res 10:2511‒2516